Structure of ternary additive hard-sphere fluid mixtures 



(N 
O 

o 

(N 

» 

o 
O 



> 

CO 

a 

C 

o 

(N 
> 
(N 

in 
in 
in 
o 

(N 
O 

c3 



13 

o 
o 



X 



Alexander MalijevskyQ and Anatol Malijevskyfj] 

Institute of Chemical Technology, 16628 Prague 6, Czech Republic 

Santos B. Yuste^j and Andres Santos^ 
Departmento de Fisica, Universidad de Extremadura, Badajoz, E-06071, Spain 

Mariano Lopez de Harcj^] 

Centro de Investigacion en Energia, UNAM, Temixco, Morelos 62580, Mexico 
(Dated: February 1, 2008) 

Monte Carlo simulations on the structural properties of ternary fluid mixtures of additive hard 
spheres are reported. The results are compared with those obtained from a recent analytical approx- 
imation [S. B. Yuste, A. Santos, and M. Lopez de Haro, J. Chem. Phys. 108, 3683 (1998)] to the 
radial distribution functions of hard-sphere mixtures and with the results derived from the solution 
of the Ornstein-Zernike integral equation with both the Martynov-Sarkisov and the Percus-Yevick 
closures. Very good agreement between the results of the first two approaches and simulation is 
observed, with a noticeable improvement over the Percus-Yevick predictions especially near contact. 
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I. INTRODUCTION 



It is widely recognized that hard-sphere fluids have 
played a key role in the development and consolidation 
of liquid state theory. For these model systems, the link 
between structural properties and thermodynamics is im- 
mediate and simple, leading to rather straightforward ex- 
pressions for the internal energy (which reduces to that 
of the ideal gas), and for the pressure equation, which 
only involves the contact values of the radial distribu- 
tion functions (rdf)|l|, 0, 0|. Nevertheless and despite the 
vast amount of literature devoted to their study, up to 
this day even the derivation of an explicit (exact) equa- 
tion of state for these systems remains as an open prob- 
lem. Under these circumstances, computer simulations 
have proved to be a useful way to derive structural and 
thermodynamic information as well as to allow the as- 
sessment of the many approximate theories proposed for 
them. These theories range from useful empirical ex- 
pressions for the contact values of the rdf or the equation 
of state to the solution of Ornstein-Zernike (OZ) integral 
equations with a given closure. And of course both theory 
and simulation increase their complexity if one considers 
mixtures rather than single component fluids, so that it is 
not surprising that the available results are much scarcer 
for hard-sphere mixtures than for pure hard-sphere flu- 
ids. In fact, only binary mixtures have received some 
attention while results for ternary hard-sphere mixtures 
and those composed of more than three components are 
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particularly limited. As far as we are aware, there is 
only one computer simulation study on the structure and 
thermodynamics of the hard-sphere additive ternary mix- 
ture Q, where diameter ratios 1:2: were considered 
at three densities and two compositions. One should 
also mention that Schaink (5) has performed a simula- 
tion study of a non-additive hard-sphere ternary mixture 
where the diameter ratios 1:1:1 were taken; the same mix- 
ture was studied by Gazzillo from an integral equa- 
tion approach. On the theoretical side, it is imperative 
to mention the pioneering work of Lebowitz j?], || , who 
solved the Percus-Yevick equation for a niulticomponent 
mixture of additive hard spheres. Also important are 
the papers by Boublfk 0, Grundke and Henderson |jTo| , 
and Lee and Levesque pdj | , in which they introduced the 
contact values, now referred to as the Boublfk-Grundkc 
Henderson-Lee-Levesquc (BGHLL) contact values, lead- 
ing to the Boubh'k-Mansoori-Carnahan-Starling-Leland 
(BMCSL) equation of state [l^] Apart from these, in 
the case of multicomponent mixtures to our knowledge 
there is only some work by Gazzillo |l3| on the thermo- 
dynamic criteria of local stability, a paper by Boublfk 
on rdf, the scaled field particle theory of isotropic hard 
particle fluids of Roscnfcld [j^ and the studies carried 
out by some of us jl6], |l^]. In these latter studies an 
interesting behavior of the rdf g%j(r) was predicted, but 
it could not be assessed in view of the then absence of 
available computer simulation data to compare with. 



On another vein, it is clear that ternary mixtures are 
typical in nature and technology. For instance, air is es- 
sentially a mixture of nitrogen, oxygen, and argon (the 
concentration of other components is much lower), and 
sea water is a mixture of H2O, Na + and Cl — . There 
is also a number of industrially important chemical re- 
actions among three components, e.g. the synthesis of 



2 



ammonia 



3H 2 + N 2 = 2NH 3 



or in ecology, e.g. 



S0 2 



1 



o 2 = so. 



Ternary mixtures of molecules whose interaction includes 
an attractive part have been studied from perturbation 
theory and van der Waals one-fluid theory [ jj), |l9|, . 

In view of the above, the aim of this paper is to pro- 
vide simulation results for hard-sphere additive ternary 
mixtures that will serve as a starting point to assess 
the accuracy and validity of some theoretical approaches. 
Specifically, we will examine five ternary systems at the 
same packing fraction and with fixed diameter ratios, so 
that they are only different because of their composition. 
Two of these cases correspond to mixtures in which the 
biggest spheres occupy over 50% of the available volume, 
followed in volume occupation by the intermediate sized 
spheres and finally by the smallest spheres. A third sys- 
tem is considered in which all species share equitatively 
the available volume, while in the last two systems it is 
the intermediate spheres which occupy the smallest vol- 
ume and either the biggest or the smallest sized ones 
follow in volume occupation. The theoretical approaches 
that we will consider will be the solution of the Ornstein- 
Zernike equation with both the Percus-Yevick [[7| and the 
Martynov-Sarkisov 21 closures and the (approximate) 
expressions for the rdf of a hard-sphere mixture derived 
in Ref. @. 

The paper is organized as follows. In order to make 
the paper self-contained, in Sect. |fj we recall the main 
results of the theoretical approaches to derive the struc- 
tural properties of hard-sphere mixtures. Section III pro- 



vides some details of the simulation and the comparison 
between simulation data and the different theoretical ap- 
proximations. We close the paper in Sect, 
discussion and some concluding remarks. 
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where hij{r) = gij{r) — 1 and Cy denote the total and 
the direct correlation functions, respectively, and jij — 
hij — dj is the series function. A general closure for the 
OZ equation may be written in the form 

dj(r) = exp [-0Uij(r) + 7„ (r) + By(r)] - 1 - 7y(r), 

(2) 

where Wy(r) is the interaction potential and is the 
bridge function. In this work we consider two approxima- 
tions to the bridge function: the classical Percus-Yevick 
(PY) theory 

B ij (r)=ln[l+ 7ij (r)]- lij (r), (3) 
and the Martynov-Sarkisov (MS) |2lJ theory 



B, 



/(r) = ^l + 2 7ii (r)-l-7ii(r). 



(4) 



We solved the OZ equation with these closures in the 
case of ternary mixtures (n — 3) using an algorithm that 
is a direct extension of the method proposed for one- 
component simple fluids |p2f . In our numerical imple- 
mentation in this paper, we used N = 2048 grid points 
with a step size Ar = 0.01. In the case of the PY closure, 
the rdf can also be obtained by numerical inversion of an- 
alytical expressions in Laplace space jl6| . Both methods 
give undistinguishable results, what gives confidence on 
the accuracy of the numerical solution of the MS closure. 

An alternative method to obtain an approximate ex- 
pression for gij (r) for a multicomponent mixture was in- 
troduced in Ref. @ . We will refer to this method as the 
RFA approach since it stemmed out of a generalization 
of a rational function approximation to structural quan- 
tities in a simple hard-sphere fluid p3[ . Working in the 
Laplace space and defining Gij(s) = J Q dr e~ sr rgij(r), 
the foregoing approach implies that G*h is assumed to 
adopt the following functional form [|l6[ glj : 



fc=i 



Ga (a) 



^^[(l+Hl-A^]- 1 ^, (5) 



II. THEORETICAL APPROXIMATIONS TO 
THE STRUCTURAL PROPERTIES OF 
MULTICOMPONENT MIXTURES OF ADDITIVE 
HARD SPHERES 



where 



Lij(s) = L^ 



L (2) s 2 



(0) 



An n-component mixture made of pi hard spheres (of 
diameter <j{) per volume unit may be characterized by 
2n — 1 parameters (for instance, the n — 1 mole fractions 
x i = Pi/p, the n — 1 size ratios Oij<J\ and the packing 
fraction r\ = ^ i=1 r?i, r\i — \pi<j\ denoting the partial 
packing fraction corresponding to species i) and involves 
n(n + l)/2 rdf gij(r). Within the usual integral equation 
approach, the OZ equation is a set of n(n+ l)/2 coupled 
equations 



p^2%k 

fc=i 



dr'h ik (\r'\)c kj (\r-r'\), (1) 



with 



Aij(s) = p., 



p=0 



-lr(2- 



p) 



ip p (x) 



E 

m=0 



(-*)' 



(7) 



(8) 



There are two basic requirements that Gij(s) must sat 
isfy. First, since gij(r) = for r < o 
(a. 



dj) /2, and the contact values 



with Oij = 
;) = finite, 



3 



this implies that (i) limg-.oo s e S(Tij Gij(s) — finite. Sec- 
ond, the isothermal compressibility kj- = finite, so that 
(ii) \im s — t o[Gij(s) — s~ 2 } = finite. Condition (i) is verified 
by construction. On the other hand, condition (ii) yields 
two linear sets of n 2 equations each, whose solution is 
straightforward: 



L 



(0) 



A + A'er, + 2X'a 



k=l 



J kj 



(9) 



L 



(i) 



A' A 
^ a ij + -^ <T i (J j + (A + X'ai)a - -a l ^ Pk^kL 

fe=i 



(2) 
kj ' 



(10) 

where A = 2tt/(1 - 77) and A' = (A/2) 2 p(cr 2 ) with (a p ) = 
^27=i Xi<T i- The parameters L^' and a are arbitrary, so 
that conditions (i) y (ii) are satisfied regardless of their 
choice. In particular, if one chooses L 



(2) 



a = 0, the 

approximation given by Eq. (|s|) coincides with the PY 
solution. If, on the other hand, we fix given values for 

9v( a tj)' we S et tne relationship l\j — 2i:aa l jg i j(a^); 
thus, only a remains to be determined. Finally, if we fix 
kt, we obtain an algebraic equation for a of degree 2n. 

In previous work with the RFA approach |ll| Q the 
BGHLL values of gijicr^) and kt given by the BMCLS 
equation of state || [l^] were considered. In this work, 
however, we will use a different approximation which was 
recently proposed by three of us |25| . Following this pro- 
posal, we assume that 



(11) 



where Zij = (aiCTj /aij)(a 2 } / (a 3 ) , and take the function 
F(r), z) to be universal in the sense that it is a common 
function for all the pairs ij. Further F is forced to comply 
with known exact relations in the point particle, equal 
size and colloidal limits. Under these circumstances, the 
simplest functional form that F may adopt is a quadratic 
function of z: 



F{rj, z) = F (r}) + F^rfiz + F 2 (r 1 )z 2 , 
where the coefficients are explicitly given by 

1 



1 — 7] 



(12) 



(13) 



F 1 (r 1 )=2(l-r,)g(a+)-^-^, (14) 

1 — 77 



F 2 (v) = ^—^-(l-2 V )g(a+). (15) 
I-77 

Here, g(u + ) denotes the contact value of the radial dis- 
tribution function of a simple hard-sphere fluid. For this 



latter, we take the one corresponding to the Carnahan- 
Starling equation of state [gfjj, namely gcs( a+ ) = (1 ~ 
?7/2)/(l - 7/) 3 . With such choice, Eqs. (0) and @ be- 
come 



1 



3 7/(1- 77/3) , 77 2 (l-77/2) 



9i ^~ l- ? 7 + 2 (I-77) 2 Zij+ (I-77) 



z 2 - 
3 y ' 



(16) 

and the compressibility factor for the mixture, from 
which kt may be readily derived, reads 

z(v) = ^bmcslW - r , (^)(^ 3 ) - (^ 2 ) 2 ) . 

(17) 

where the compressibility factor associated with the BM- 
CSL equation of state [g, |l^| is given in the present no- 
tation by 



(1-77)2 ((T 3)2 



^BMCSL 



(v) 



1 



St,^}^ 2 ) , 77 2 (3 - 77)(a 2 } 



2\3 



1-r] {l-7]) 2 (a 3 ) (1 



r]) 3 {a 3 ) 2 



(18) 



Equation jl^ ) represents in general a significant improve- 
ment over the BGHLL contact values |EaL On the 
other hand, the BMCSL equation of state (|18[) performs 
slightly better than Eq. (|17|). Although the RFA can be 
implemented by making any choice for gij{crfj) and kt, 
here we have taken, in addition to the contact values ( jig) , 
the isothermal compressibility associated with Eq. ( |l7j ) 
in order to enforce thermodynamic consistency. 



III. SIMULATION DETAILS AND RESULTS 

We used the standard NVT-Monte Carlo method with 
periodic boundary conditions, employing a cell index al- 
gorithm with six different cell sizes corresponding to a 
number of interactions. The simulation cubic box con- 
tained N — 2700 particles in each case but one (case C), 
where N — 6777 particles were used. 

The initial system with no overlaps was generated by 
random insertion of particles to an originally empty box. 
The sequence we used is the following: the largest parti- 
cles were inserted first and the smallest ones at the end. 
Particles were mixed during this procedure. Starting 
with this initial configuration, we generated the Monte 
Carlo chain as follows. The acceptance ratio of trial 
moves was adjusted to 10-15% for all the components. 
Each run was divided into 21 blocks, each of which in- 
cluded about 10 9 of the equilibrium configurations gen- 
erated and contained 300-500 analyses of the calcula- 
tion of the rdf gij(r) in the whole range of 1200 intervals 
ri ± Ar/2 (where the step size was Ar = 5 x 10~ 3 (7i) up 
to a distance 601. The analysis was performed after 1000 
trial moves per particle (more precisely after 1000A r trial 
moves of a randomly chosen particle). The first block was 
then discarded and the next 20 were used to sample the 
configuration space, calculate mean values for the entire 
run and estimate the errors. 



The systems we examined had the same packing frac- 
tion r\ = 0.49 and fixed diameter ratios O2/01 = 2 and 
03 1 <J\ = 3 (for convenience and without loss of generality 
we have chosen the value of the diameter of the smallest 
spheres to be always 1) so that their only difference lies 
in the composition. They are identified as 

(A) x x = 0.7, x 2 = 0.2, x 3 = 0.1, 

th/ti = 0.14, m/v = 0-32, 773/77 = 0.54, 

(B) xi = 0.6, x 2 = 0.2, x 3 = 0.2, 

7/1/77 ~ 0.08, % /7/ ~ 0.21, 773/77 ~ 0.71, 

(C) X\ — , X 2 — 25J, Xz = 2§i, 

m/v = m/v = m/v = |, 

(D) xx = 0.85, x 2 = 0.05, x 3 = 0.10, 

771/77 ~ 0.22, 7/2/7/ ~ 0.10, 773/77 ~ 0.68, 




,0.0 



<2> 



C5- 



to 



V 



(E) xx = 0.90, x 2 = 0.07, x 3 = 0.03, 
7/1/77 ~ 0.396, 772/77 ~ 0.247, 7/3/7/ 



0.357. 



These systems have been located in two different 
ternary diagrams, one with respect to mole fractions and 
the other one corresponding to partial packing fractions, 
shown in Figs. ^ and |^, respectively. In these diagrams 
we have also included the two systems with diameter 
ratios 02/eri = 2 and <j 3 /(Jx — 10/3 that were studied 
by Sindelka and Boublfk Q and which we have labeled 
SB1 (xx = x 2 = x 3 = \; 771/7/ ~ 0.022, 7/2/77 ~ 0.174, 
773/77 ~ 0.804) and SB2 (xx = \, x 2 = |, x 3 = \; 
771/77 ~ 0.054, 772/7/ ~ 0.285, 773/7/ ~ 0.661). It should be 
pointed out that cases A and B (as well as the systems 
SB1 and SB2) correspond to the situation 7/1 < 7/2 < 7/3, 
while in case D one has 772 < 7/1 < 7/3, in case E 
V2 < Tj3 < 771 and in case C 771 =772 =773. This, in 
our view, allows us to examine the very different situa- 
tions that arise depending on which species occupies the 
largest volume. 

The results of our calculations, both theoretical and 
from the simulations, are displayed in Figs. ||-|8|. In Fig. 
U we show the contact values for all five systems as func- 



tions of the parameter Zij defined below Eq. (|1 1|) . In 
this instance we have considered the PY, BGHLL, MS, 
and scaled particle theory (SPT) values (ll| |7j, as well 
as those given by Eq. ([l6|). In the case of the MS ap- 
proximation we actually get a set of points that have 
been joined by a line interrupted at Z33 = 1.481 (case 
D) since there is no convergence in cases C and E. The 
fact that this line is a smooth one shows that the numer- 
ical values obtained from the MS approximation seem to 
be consistent with the "universality" assumption (11). 
The comparison with the simulation data indicates that 
for z > 1 the new proposal, Eq. (|l6|), improves over the 
BGHLL prescription (while for z < 1 it is only slightly 
worse) and both are clearly superior to the SPT recipe. 
The PY values are very poor, while the MS approxima- 
tion tends to underestimate the contact values for z > 1. 
This provides some support to the use of Eq. (Jlq) and 



FIG. 1: Ternary diagram showing the mole fractions of the 
five cases A-E considered in this paper, as well as the two 
cases (SB1 and SB2) considered by Sindelka and Boublfk @. 
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FIG. 2: Ternary diagram showing the (relative) packing frac- 
tions of the five cases A-E considered in this paper, as well 
as the two cases (SB1 and SB2) considered by Sindelka and 
Boublfk |). 



Eq. ( |l7j ) (this latter to compute Kt) within the approxi- 
mate scheme to derive the rdf Oy (r) for ternary mixtures 
that was introduced in Ref. |16 and briefly sketched in 
the previous section. 

Figures |||-|| show all the rdf gij(r) = 1,2,3) as 
functions of the shifted distances r — Uij for the five differ- 
ent systems considered. Also included in these figures are 
insets with an enlarged scale around gij (r) = 1 in which 
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FIG. 3: Plot of the contact value gij(cr^) as a function of the 
parameter Zij = (aiaj)(a 2 ) / (a 3 ) for ternary additive hard- 
sphere mixtures at a packing fraction r\ = 0.49. The circles 
are simulation data for the five cases A-E considered in this 
paper. The lines are theoretical predictions: Eq. ([L6|) ( — ), 
BGHLL ( ), SPT ( ), MS ( ), and PY (■ ■ ■ )■ 



we have plotted gij(r) versus the actual rather than the 
shifted distance. It should be noted that, as said before, 
the solution to the OZ equation with the MS closure did 
not converge for cases C and E. This is a consequence 
of the fact that a term under the square root in Eq. (^) 
becomes negative at high densities; the authors of the 
MS closure speculate that the lack of convergence may 
be a signal of phase transition S| . The analysis of this 
conjecture is beyond the scope of this paper. The study 
of Figs. ^-|| indicates the following. The RFA approach 
provides an excellent overall agreement with the simula- 
tion results, which is especially good in the region around 
contact. Something similar occurs with the solution to 
the OZ equation with the MS closure, except that this 
solution tends to underestimate the contact value of 523 
and 333. The PY closure clearly yields the poorest results 
especially in the region around contact. All three theo- 
retical approaches lead to almost identical results beyond 
the first minimum and exhibit a rich fine structure as was 
already pointed out for another ternary system in Ref. 
]l6| . The fact that the simulation results also exhibit this 
structure is in our view remarkable. It should be noted 
that there are slight quantitative differences around the 
first minimum, which is more pronounced in the theoret- 
ical solutions than in the simulation. Except for cases C 
and E where the fine structure is rather similar, in the 
other three cases the fine structure is case dependent. As 
may be observed in Fig. 0, cases C and E correspond to 



a situation where all partial packing fractions are rather 
similar. Interestingly enough, when this happens, i.e. no 
species is dominant with respect to volume occupation, 
all the rdf gij (r) relax to 1 following an ordered sequence 
of damped oscillations. Finally, it is also worth mention- 
ing that, for a given system, the form of the fine structure 
of the gij{r) is almost the same for all pairs. In fact, such 
fine structure seems to evolve smoothly as <Jij increases 

(an = 1, <7i 2 = 1-5, (7 2 2 = CTi3 = 2, (723 = 2.5, (7 33 = 3) 

as one can easily see by following the sequence top pan- 
ncl left , top pannel right, middle panncl right, middle 
panncl left, bottom panncl right, bottom panncl left in 
Figs. |, |, and 0. 



IV. CONCLUDING REMARKS 

The results of the previous section deserve some further 
comments. As far as the simulations are concerned, we 
have here provided further data on the thermodynamic 
and structural properties of ternary additive hard-sphere 
fluid mixtures that extend and complement those of Ref. 
Q. The quality and reliability of these new data is re- 
flected in their ability to capture the rich fine structure 
that had been observed earlier in connection with the 
RFA approach |[(|. Both the RFA results and the ones 
derived from the solution of the OZ equation with the 
MS closure are in very good agreement with the simula- 
tion data, but the latter give less accurate contact val- 
ues, the OZ equation needs to be solved numerically, and 
it presents convergence problems when the partial pack- 
ing fractions of all three species have similar values. In 
any case these two theoretical approaches do represent 
an improvement over the Percus-Yevick theory. Finally, 
we have only carried out a preliminary qualitative anal- 
ysis of the rich fine structure that arose in the systems 
we examined. By restricting ourselves to a fixed total 
packing fraction and given diameter ratios, we attempted 
to investigate the effect of partial volume occupation by 
each species on the resulting structure. It thus appears 
interesting to assess the effect of different total packing 
fractions and (or) size ratios. We may address these and 
other related issues in multicomponent systems in the 
future. 
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FIG. 4: Radial distribution functions gij(r) for a ternary mixture with diamters <7i = 1, 01 = 2, 1T3 = 3 at a packing fraction 
r\ = 0.49 in case A [x\ = 0.7, £2 = 0.2, x-$ — 0.1). The circles are simulation results, the solid lines are the RFA predictions, 
the dotted lines are the PY predictions, and the dashed lines are the MS predictions. 
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FIG. 5: The same as in Fig. kl but in 
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B (xi = 0.6, x 2 = 0.2, x 3 = 0.2). 
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